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Abstract 

We numerically study a disordered version of the model for DNA 
denaturation transition (DSAW-DNA) consisting of two interacting 
SAWs in 3d, which undergoes a first order transition in the homoge- 
neous case. The two possible values eUr and ice of the interactions 
between base pairs are taken as quenched random variables distributed 
with equal probability along the chain. We measure quantities aver- 
aged over disorder such as the energy density, the specific heat and 
the probability distribution of the loop lengths. When applying the 
scaling laws used in the homogeneous case we find that the transi- 
tion seems to be smoother in the presence of disorder, in agreement 
with general theoretical arguments, although we can not rule out the 
possibility of a first order transition. 

PACS: 87.14Gg,87.15.Aa,64.60.Fr,82.39.Pj 

1 Introduction 

The DNA denaturation transition is an old-standing open problem [1] and 
one finds in the literature a large number of different models that look in de- 
tail at various aspects [2]-[5], for instance reproducing the double-helix struc- 
ture and the corresponding denaturation behaviour in the whole temperature- 
torsion plane [I] or taking into account to different extents the role of the 
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stacking energies [31 IS]- From the experimental point of view, one observes a 
multi-step behaviour in light absorption as a function of temperature which 
suggests a sudden sharp opening of clusters of base pairs in cooperatively 
melting regions. Therefore one expects that a theoretical model, correctly 
reproducing the experimental behaviour, should undergo a sharp transition 
and there has recently been a lot of interest in models possibly displaying a 
first order transition in the homogeneous case [l]-{llj. 

Most programs, for example MELTSIM [12], predict the experimentally 
observed melting curves using a model originally introduced by Poland and 
Scheraga [21 [11] that takes into account the different entropic weights of 
opened loops and double stranded segments. In the first studies of this model, 
excluded volume effects were completely neglected and a smooth second order 
transition was predicted in two and three dimensions. By solving the model 
including the entropic weights of self-avoiding loops [13], one finds a sharper 
but still second order transition. In this way the self-avoidance between 
bases within the same loop is taken into account, but other mutual excluded 
volume effects are still neglected. 

In a previous work [7] we pointed out the importance of excluded volume 
effects between different segments and loops by introducing and numerically 
studying a model inspired by the former of Poland and Scheraga, consisting 
of two interacting self-avoiding walks on a 3c? lattice (SAW-DNA). 

In the limit of infinite chain length, the SAW-DNA model undergoes a first 
order phase transition from the double strand to the molten single-stranded 
chain state when varying the temperature. The order parameter, which is 
the energy or the density of binded base pairs, varies abruptly from one to 
zero. 

It has been theoretically demonstrated that the transition in the homoge- 
neous model by Poland and Scheraga is of first order when excluded volume 
effects are completely taken into account, using the entropic weight of a self- 
avoiding loop embedded in a self-avoiding chain [HI El- This can be obtained 
as a particular case of a more general result on polymer networks [15] . It has 
also been analytically shown that only considering the self-avoidance between 
the two different chains leads to a first order transition [TU1 [TO] . 

Moreover, several numerical investigations [H [161 El EEH] of the homoge- 
neous SAW-DNA model carefully measured the probability distribution of 
the lengths of denatured loops P(l) at the critical point in 3d and 2d. These 
confirmed the theoretically expected power law P(l) oc l/l c with exponent 
c > 2 (in agreement with the transition being first order) in both dimensions. 
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In the same SAW-DNA model, a first order unzipping transition is pre- 
dicted in the presence of an external force that pulls apart the two DNA 
strands [H]. The behaviour in the temperature-force plane is numerically in- 
vestigated in [19], and a phase diagram with a re-entrance region is observed. 

The version of the Poland-Scheraga model on which programs for pre- 
dicting experimental denaturation curves are based, contains another main 
parameter, the cooperativity factor a, which takes into account the activa- 
tion barrier to open a loop and has the effect of sharpening the transition 
when considering finite size chains. It was shown [20] that by using the value 
c ~ 2.15 that characterizes P(l) in 3d (instead of the usual value c ~ 1.75, i.e. 
the exponent of an isolated self-avoiding loop), and by slightly correspond- 
ingly varying the cooperativity factor, one still obtains a multi-step melting 
behaviour well in agreement with experimental results. Nevertheless, the 
relevance of self-avoidance for the experimental DNA denaturation is still an 
open question [21], [22] . 

In this work, we present numerical results on the denaturation transition 
in the interacting SAW model in 3d in the presence of quenched disorder 
(DSAW-DNA). The main aim of our study is to attempt to clarify whether 
disorder is relevant, therefore we introduce a model as simple as possible, 
since we expect that slightly more realistic versions should behave similarly. 
Studies in the literature on the effects of sequence heterogeneities on other 
simple models for DNA denaturation [2U] 12"^] neglect self- avoidance. Apart 
from its importance to experimental melting, we find that this is an intrigu- 
ing statistical mechanics model in itself. On the one hand, from general 
theoretical arguments [25]- [27], one may expect that the transition should 
become smoother in the presence of disorder, although this is a peculiar kind 
of first order phase transition, corresponding to a tri-critical point in the 
fugacity-temperature plane, and it is characterized both by an a = 1 specific 
heat exponent and by a diverging correlation length [7]. On the other hand, 
numerical results on P(l) at the critical temperature for a single disordered 
sequence in this model [9] seem to show that the order of the transition does 
not change, as also suggested from the topological considerations explaining 
the sharp transition in the homogeneous case. 

We shall consider the simple case in which there are only two possible 
base pair interactions, i.e. €at which describes the Adenine-Thymine couple 
and ice for the Guanine- Cytosine one, with ice = 2eUr- To be realistic, one 
should choose closer values but, if disorder is relevant, they could make it 
difficult to find numerical evidence for its effect, being reasonable to expect 
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that the disorder possibly changes the order of the phase transition as soon 
as £at 7^ £gc- Therefore we assumed to take ice — ^at as an interesting 
compromise. 

Moreover, we take the interactions to be independent quenched random 
variables, identically distributed with equal probability 1/2. Again, it should 
be noted that in more realistic models the interaction energies are chosen to 
depend also on the first neighbours (to take into account at least partially 
the stacking energies) and that there are probably long-range correlations in 
the base pair distribution, the ratio of the GC to AT content being in any 
event a highly varying quantity usually far from 1. We are neglecting both 
of these effects, but we find that our simplified model does already display 
an intriguingly rich behaviour. It seems to be a useful starting point for 
understanding the thermodynamical properties of this kind of systems in the 
presence of disorder and their relevance for describing experimental DNA 
denaturation transitions. 

2 Model and Observables 



Let us define two A-step chains with the same origin on the 3-d lattice, uj 1 = 
{c^o, . . . , u^} and u 2 = {uq, . . . , uj%}, with ouf G Z 3 and u\ = Uq = (0, 0, 0). 
The Boltzmann weight of a configuration (a; 1 ,^; 2 ) of our system is 




(1) 

where the {e{\ are independent quenched random variables, identically dis- 
tributed according to the bimodal probability 

V(e) = \ [S(e - e AT ) + 8(e - e GC )\ • (2) 

Thermodynamic properties of the system only depend on the reduced 
variables eAr/ksT and 1Z = £at/^gc: the homogeneous case corresponding 
to 1Z = 1. Here we will take Iat = 1 and ice — 2, i.e. 1Z = 1/2. For the 
sake of clarity, we also fix the Boltzmann constant to = 1, measuring the 
temperature in £at units. 

Let us introduce the different observables in the well studied homogeneous 
case. In the thermodynamical limit, the transition is a tri-critical point in the 
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fugacity-temperature (z — T) plane [7], described by the crossover exponent 
0. One has: 

z - z c oc (T - T c f (3) 

The order parameter characterizing the transition is the density of closed 
base pairs, that behaves like the energy density e = E N /N, where [7J 125] : 

f 1/(T-T C ) T^T C + 
E N (T) ~ J iV* T = T C (4) 

[ iV^-T) 1 ^- 1 T^T~ 

Therefore the value = 1 of the crossover exponent corresponds to a first 
order transition, in which e goes discontinuously (in the thermodynamical 
limit) from the zero value of the high-temperature coiled phase to a finite 
value at T c . We stress again that it is a peculiar kind of first order transition 
with a diverging correlation length and absence of surface tension (the prob- 
ability distribution of the energy is nearly flat at the critical temperature). 

In the homogeneous case, both when self-avoidance is completely taken 
into account (i.e., a first order transition with = 1 in d = 3 and d = 2) 
and when it is neglected (the random walk model, which undergoes a second 
order transition with 0=1/2 in d = 3 and a first order transition for d > 5), 
one finds that the finite size behaviour of different quantities is in agreement 
with given scaling laws. The total energy En, its probability distribution 
P(En), and the maximum of the specific heat c™ ax behave as [3 [28] : 



E N (T)/N* = ~h[(T-T c )N*] (5) 
c max (N) oc N 2 ^ 1 (6) 
P N (E)N* = f(E/N+) at T = T C , (7) 

where / and h are scaling functions. 

Moreover, one can get an independent evaluation of the crossover expo- 
nent from the probability distribution of the loop lengths, that at T c is in 
agreement with the law: 

P(l) oc p (8) 

where the exponent is related to <fi through the equation <fi = min{l, c — 1}. 
In the 3d SAW-DNA, one finds numerically [TU [IB] c = 2.14(4), in perfect 
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agreement with the theoretical prediction [5] c ~ 2.115, that implies a first 
order transition. 

It is also interesting to note that P(l,T) is more generally expected to 
behave as [18]: 

and that the correlation length £(T) should diverge when approaching the 
transition point with the power law: 

( \T — T c \~ l for c > 2 

£(T) oc | |T _ £|-i/ (c -i) for 1 < c < 2. (10) 

In the presence of the independent quenched random variables e = {e^}, 
one has to introduce quantities averaged over disorder: 



where 



E N (T) 


= E N , e (T) 


(11) 


CN,T 


= C7V, e (T) 


(12) 


Pn{E) 


= PnAE) 


(13) 


Pn(1,T) 


= PnA^T), 


(14) 


0(T) = J 


r deV(e)O e (T), 


(15) 



the numerical evaluation of these quantities being as usually performed by 
averaging over a (large) number of different disordered configurations (sam- 
ples). One should note that the energy density e = — uat^at — u gc^gc 
is now quantitatively different from (minus) the number of contact density 
n = riAT + n GCi though it seems reasonable to expect that the two quantities 
behave similarly. 

In particular, we will look at the previously discussed scaling laws 
on the averaged energy, the averaged specific heat maximum, the averaged 
probability distribution of the energy and the averaged probability distribu- 
tion of the loop lengths. 



3 Simulations 

We used the pruned-enriched Rosenbluth method (PERM) [29] , with Marko- 
vian anticipation [30], which is particularly effective to simulate interacting 
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polymers [21] • Moreover, we used an ad hoc bias for the present model. 
When the second chain has to perform a growth step and the end of the first 
one is in a neighboring site, instead of doing a blind step multiplying the 
weight by a factor e~ ei l kBT if the new contact is formed, we favor contacts 
choosing the step towards the end of the first chain with an appropriately 
higher probability and correcting the weight accordingly 

For each considered sequence we have performed 16 independent runs 
at different temperature values with 3^-5 million independent starts for 
each run. This turns out to generate enough statistics up to chain lengths 
N = 800, since the number of times that the system reaches the largest N 
value in the simulation is of the order of the number of independent starts. 

We compute the behaviour of the energy, of its probability distribution 
and of its derivative (i.e., the specific heat) as a function of temperature by 
re- weighing the data at the chosen temperatures of the set. The statistical 
errors on the values for a given sequence are evaluated using the Jack-knife 
method and we checked that they are definitely smaller than the errors due 
to sample-to-sample fluctuations, which are our estimate of the errors on 
averaged quantities. 

Moreover, we checked thermalization and the correct evaluation of the 
errors, particularly in the case of P(l,T), by comparing results obtained at 
two different sets of temperatures (see the Table) with two different meth- 
ods. In the first case two copies of the systems evolved simultaneously and 
independently (the algorithm being accordingly modified) and the computed 
quantities are practically evaluated from two independent runs of 3 millions 
of starts each. In the second set one copy performed 5 millions of indepen- 
dent starts. We obtained perfectly compatible results from the two sets of 
simulations. 

We considered 128 different samples. We note that the disorder configu- 
rations for different N values are not completely independent, nevertheless it 
seems reasonable to expect to observe the same scaling properties. Therefore 
we find our statistics to be sufficient for giving a first insight onto the be- 
haviour of various quantities and we also stress that the crossover exponent 
can in principle be obtained from data on Pn(1) at T c corresponding to a 
single N value. The whole simulation took about 15000 hours (on COMPAQ 
SC270 and, to a lesser extent, on Forshungszentrum Jiilich CRAY-T3E). 
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rpa 

i 


0.8 


0.95 


1.08 


1.125 


1.15 


1.175 


1.2 


1.3 


rpb 

% 


0.875 


1.015 


1.05 


1.1 


1.13 


1.16 


1.19 


1.25 



Table 1: The considered temperatures for the two sets of simulations per- 
formed for each sample. In the first case (a) two copies of the system evolved 
independently and statistics were collected over 3 millions of starts, whereas 
in the second case (b) only one copy was simulated and statistics were col- 
lected over 5 millions of starts. 

4 Results 

We shall focus on results and compare them to what is observed in the ordered 
case. 

First of all we find, at least for some sequences (about one third of the 
sequences for the largest length), the expected multi-step behaviour of the 
energy density (see Fig. [TJ and correspondingly several peaks of the specific 
heat (see Fig. that behaves qualitatively as the derivative of the density of 
closed base pairs with respect to the temperature, i.e. the differential melting 
curve which is usually experimentally measured. 

Qualitatively speaking, the presence of multi-steps becomes more prob- 
able with increasing N, whereas the temperature region in which this be- 
haviour is present becomes narrower. This is in agreement with the ob- 
servation that at the critical point the usual argument for proving the self- 
averaging property of the densities of extensive quantities, such as the energy 
density, fails since one can not divide the system in nearly independent sub- 
systems due to the diverging correlation length. This is true also in the 
SAW-DNA where, despite the first order transition of the homogeneous case, 
one finds a diverging correlation length. On the considered N range we ob- 
serve a strong sample- dependent behaviour, as in the case also of -P/v,e(0 at 
T ~ T c . It is therefore possible that the results of the application of the 
scaling law of the homogeneous case to averaged quantities have to be taken 
with some care [33] . 

In Fig. [3] we plot our data on the disorder averaged specific heat for the 
considered chain lengths. This shows that, also in the presence of disorder, 
the maximum of the specific heat increases as a function of the chain length. 
Nevertheless, as discussed in detail in the following, it seems to behave as 
jV 2 ^ -1 with exponent smaller than one. 
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Figure 1: The energy density e^^iT) for two particular sequences with 
N = 800, one of which displays a two-step behaviour. Note the plateau at 
^N,e(T) = in the coiled phase. The plotted region is the one around T c 
in which the energy varies more rapidly and there are evident differences in 
the behaviour from sample to sample. At lower temperatures quantities are 
nearly sample-independent and the energy density slowly decreases towards 
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Figure 2: The specific heat Cw )e (T) for two different disordered sequences 
with N = 800 (the same sequences as in previous figure). Also in this case 
there is qualitative agreement with the experimentally observed behaviour 
in the differential melting curves, though the model should be improved in 
order of really comparing. For instance, particularly for short sequences, the 
plateau one gets in the coiled phase is usually higher than the low tempera- 
ture one because of residual stacking energies which is an effect completely 
neglected here. 
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Figure 3: The disorder averaged specific heat for the considered chain 
lengths. 



11 




-1.4 1 1 1 1 1 1 1 1 1 1 1 

0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 1.2 1.25 1.3 



T 

Figure 4: The disorder averaged energy density, for the considered chain 
lengths. 
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Then in Fig. H] we present our data on the disorder averaged energy den- 
sity. We note that curves corresponding to different N values cross roughly 
at the same temperature T c ~ 1.15 within the errors (see the insert), there- 
fore suggesting a transition still of first order. Actually this should mean 
a jump of the energy density in the thermodynamical limit, from zero for 
T > T c to the finite value e(T c ) ~ —0.5. Nevertheless, as displayed in Fig. 
[5^, the expected scaling law Eq. (jSJ), which is well verified in the homoge- 
neous case at least in the region T > T c , is here definitely not fulfilled with 
0=1. When fitting the data according to this scaling law (see Fig. [5b) one 
finds a slightly higher critical temperature value T c = 1.155 and a crossover 
exponent = 0.8475. 

Since we were looking for a second order transition, i.e. for an exponent 
< 1, we required in the fit that data obey the scaling law on both sides 
of the critical temperature (for the random walk model in 3c? [7], where the 
transition is of second order, a good scaling was observed on both sides of 
T c ) . In any event we stress that by fitting data only in the high-temperature 
region T > T c one would get a definitely lower value of the crossover exponent 
~ 0.6. 

We note that the data roughly agree with the law but there are still 
corrections to scaling. We also checked that the disorder averaged number of 
contact density displays a similar behaviour (i.e. it follows the same scaling 
law). 

In Fig. [6] we present our data on the specific heat maximum, that would 
agree well with an exponent ~ 0.85 (our fit to the expected behaviour Eq. 
(ED gives = 0.815). However the data would also be compatible (within the 
errors) with = 1, if we were to neglect the smallest chain length N = 100. 

We tried, following [33J, to analyze the data on the energy densities and 
the specific heats in terms of a reduced sample dependent temperature T — 
T*(N), where T*(N) is defined as the temperature for which the specific heat 
for sample i and chain length N reaches its maximum. This analysis leads 
to the same conclusion as the conventional one, that the disorder averaged 
energy density does not scale with a = 1 crossover exponent and the 
disorder averaged specific heat maximum seems to diverge more slowly than 
linearly with N when taking into account all the considered chain lengths. 
We computed the average distance between the sample dependent critical 
temperatures dT c {N) = \ T c( N ) ~ Tj(N)\ (W = 128 being the 

number of samples) which seems to go to zero for increasing chain lengths 
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Figure 5: On the top we plot the disorder averaged energy density for the 
considered chain lengths as a function of N(T — T c ) with T c = 1.15. On the 
bottom we present the scaling law with <fi = 0.8475 and T c = 1.155. 
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more slowly than 1/N, again suggesting a second order transition. 
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N 

Figure 6: The behaviour of the specific heat maximum compared with a fit 
to the exponent = 0.815 and to = 1. 



We looked at the probability distribution of the energy density at the 
critical temperature T — 1.15 ~ T c (see Fig. Uh)- In the homogeneous case 
it displays a large and flat scaling region extending to a value e* which do 
not depend on N, with deviations from the scaling law Eq. (j7|) for e < e* 
[7J. Here we find that the scaling law with = 1 is not well fulfilled. For 
the sake of comparison we also plot (in Fig. [7b) the behaviour of our data at 
the slightly higher temperature value T = 1.155 with a crossover exponent 
= 0.8475. 

Let us now consider data on the disorder averaged probability distribution 
of the loop lengths Pn{1) which are plotted in Fig. [8] at T = 1.15 ~ T c . 
We checked that the behaviour is very similar when temperature slightly 
varies. Again, a fit to the expected power law oc l/l c of the whole data set 
suggests an exponent c < 2, i.e. a second order transition. Nevertheless, a 
closer inspection of the behaviour makes evident some curvature and when 
restricting the fit to the / range 1 << I << N (which is also the region in 
that the hypothesis Pn(1) oc 1// c should be better verified) we find larger c 
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Figure 7: On the top we plot the disorder averaged probability distribution 
of the energy density at T = 1.15 ~ T c , Pjv(e) as a function of — e, which 
corresponds to the particular case of the corresponding scaling law with = 
1, since P/v( e ) = N Pn(E/N). On the bottom we present the scaling law, 
P N (E)N^ as a function of -E/N*, at T c = 1.155 with = 0.8475. 
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values and the data for N = 100 alone are definitely consistent with a first 
order transition. 
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Figure 8: The disorder averaged probability distribution of the loop lengths 
at T = 1.15 ~ T c for N = 100 and N = 800, compared to the behaviour a/l c 
with c = 2. 



To perform a more quantitative analysis, in Fig. [9] we consider [34] the 
momenta < l p >, which are expected to behave as N Up . In particular the 
combination In < l p >/ In N is expected to be linear in p on a large range of 
p values and one should be able to evaluate c — 1 by extrapolating the linear 
behaviour down to p = 0. We effectively observe a quite linear behaviour for 
p > 1 but by extrapolating we get the definitely different values c = 2.07 for 
N = 100 and c = 1.91 for N = 800. On the one hand data for the longest 
chain length should be the most meaningful, therefore confirming a second 
order transition; on the other hand it should be stressed that P/v j(E (7) is a 
difficult quantity to measure and we can not rule out the possibility that our 
statistics are inadequate for correctly evaluating it for the largest considered 
chain lengths. 

Nevertheless, the observed curvature on the log P^(l,T c ) vs\og(l) be- 
haviour and the corresponding strong ^-dependence in the obtained esti- 
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Figure 9: Data on a p = ln(< l p >)/lnN, from data on P(l) at T = 1.15 ~ 
T c , plotted as function of p. The two fits give c = 2.07 and c = 1.91 respec- 
tively. We note that the linear behaviour is no more satisfied for p > 6 in 
data corresponding to the largest iV values. 
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mations could be related to finite size effects. One should note in particular 
that in the case of a second order transition the finite size corrections to the 
critical temperature (i.e. T C (N) — T c (oo)) are expected to approach zero in 
the thermodynamical limit more slowly than 1/N. Therefore, it could be 
misleading to fit the data according to the power law P(l) oc l// c which is 
expected to be fulfilled only exactly at T c . Moreover we have a large un- 
certainty on the evaluation of T c itself. For these reasons, in order to gain 
further insight into the behaviour of the system, we also attempt to compare 
the measures at the different temperatures considered with the more general 
law Pjy(l,T) oc exp(— 1/£n(T))/I c . This should allow to take partially into 
account finite size effects, by introducing a possible finite correlations length 
£jv(T) which measures both the distance from the effective T C (N) and the 
finite size corrections to the thermodynamical limit correlation length. 

Interestingly enough, in the temperature region T > T c ~ 1.15 one gets 
1/£at(T) compatible with zero within the errors (it can be also negative, 
though small) apart from the shortest chain lengths. In particular, data for 
N = 100 are consistent with a detectably finite correlation length for all the 
temperatures studied. The corresponding evaluation of the exponent c is in 
this case smaller than that found when the effect of the finite correlation 
length is neglected, giving usually c < 2. 

In detail, we performed three-parameter fits of data on P/v(/,T) at dif- 
ferent temperatures by neglecting both the very first and the last I values. 
In particular the dependence on the number of neglected initial points has 
been studied by restricting the range to I > 5 and to I > 2. Whenever the 
fit gives a l/£jv(T) < unphysical value, the corresponding c is estimated 
from the power law behaviour, i.e. by imposing 1/£at(T) = 0. Moreover we 
also considered the fits to the power law in all the cases in which an inverse 
correlation length compatible with zero turns out, again by looking only at 
the / > 5 range or at the whole I > 2 one. Finally, we performed linear 
extrapolations to the 1/N — > limit. 

Despite the introduction of the correlation length, the evaluations of c 
seem generally to depend on temperature and on the number of initial points 
included in the fit, and its extrapolations turn out to be compatible with a 
first order transition. Nevertheless, the values obtained in the region T ~ T c , 
where the largest sizes can be fitted to the power law and which are there- 
fore expected to be the most significant results, would definitively suggest 
a second order transition, characterized by a quite large crossover exponent 



19 



~ 0.9. It should be noted that in this region one gets well compatible 
estimations from the different methods we used. In conclusion, the analysis 
is consistent with the initial qualitative observation on the P(l) behaviour 
at T ~ T c , that the transition seems to be of second order when considering 
the longest chain lengths most meaningful. 

It would be interesting to perform on P/v(Z, T) an analysis by introducing 
rescaled temperature but we unfortunately are not able to re-weigh data at 
different T values in this case. To study whether this would not affect the 
qualitative behaviour, thereby confirming a second order transition, is left 
for future investigation. 
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Figure 10: The results on the inverse correlation length 1/£(T) behaviour of 
the fits of data on P/v(7, T) as obtained by different methods: a) For each N 
and temperature the values are obtained from a three-parameter fit to the 
law oc exp(—l/^(T))/l c neglecting only the very first (in particular the first 
two) and the last / values, b) The same as in a) but restricting the considered 
range to I > 5. In both cases values for different N are linearly extrapolated 
to 1/N — *■ in order to give the plotted estimations. We stress that because 
of the fits involved the given errors are to be considered only very indicative, 
our uncertainty on the estimations being better expressed by the fluctuations 
between values obtained with the different methods. 
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Finally, we present in Fig. [TU] data on the inverse correlation length, 
1/£(T) as evaluated by extrapolating linearly to the limit 1/N — > values 
for the different chain lengths (here we only consider the temperature range 
where it is definitely larger than zero). The obtained qualitative behaviour 
suggests again a second order transition, since 1/£(T) seems to go to zero 
more rapidly than (T c — T). A fit to the expected power law behaviour 
l/f(T) oc (T c - T) 1 /^- 1 ) would give c ~ 1.7. We checked (on the N = 800 
data) that the qualitative behaviour of £(T) does not change if we were to 
fit the data on P/v(/,T) in the low temperature region according to the law 
<Q by imposing c = 2 (or c = 2.1). 

As a last remark, we stress that our statistics are unfortunately inad- 
equate for performing a more quantitative analysis, particularly on £(T). 
Actually, because of the algorithm we are using, data on P(l,T) in the low 
temperature range are reliable up to smaller I values rather than in the re- 
gion T > T c (i.e. very long bubbles are usually not well sampled at low 
temperatures due to their negligible probability). Once again, we are led to 
the conclusion that the most significant results on P(l, T) should be those 
obtained in the T ~ T c region for the largest chain lengths, which is our best 
numerical evidence for the transition being of second order. 

We have some preliminary numerical results on our model with different 
values of TZ, where the crossing of the energy densities becomes less evident 
the more different the energies e^r and €gc are - The transition seems to 
be characterized by a varying crossover exponent which becomes smaller as 
1Z diminishes. In particular, by applying the scaling law to the disorder 
averaged energy densities we get ~ 0.8 from data for €at = 1 and ecc = 4, 
and ~ 0.7 from data for 6at = and ecc = 1- The P(l) behaviour 
at the critical temperature would suggest similar values too, though also in 
these cases one gets different results by restricting the range to the region 
1 << / << N. Here we only would like to mention that an alternative 
explanation for these results is that one is always looking at a first order 
transition but that the scaling laws of the homogeneous case are not well 
fulfilled anymore [33]. More extensive simulations would be necessary in 
order to clarify this issue. 
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5 Conclusions 



Summarizing, we studied a disordered model for DNA denaturation transi- 
tion consisting of two interacting self-avoiding chains in which we chose the 
two possible values of the interaction energy to be ecc = ^at, distributed 
with equal probability. Despite of the extensive numerical simulations per- 
formed, it is difficult to definitively discriminate between a first order (as in 
the homogeneous case) and a smoother second order transition. 

As a matter of fact, we find that the energy densities as a function of 
temperature for the different chain lengths considered roughly cross within 
the errors therefore suggesting a discontinuity of this quantity in the ther- 
modynamical limit. However, the application of the scaling laws that are 
verified in the ordered case indicate a smoother second order transition with 
strong corrections to scaling. 

When looking at the averaged probability distribution of the loop lengths, 
one gets again a crossover exponent </> definitely smaller than one from the 
data corresponding to the largest considered chain lengths. Nevertheless, 
data display some curvature also at T ~ T c and higher values of the exponent 
are obtained when restricting the range to 1 << I « N. A momenta 
analysis of this probability distribution gives c < 2 (i.e. <p < 1) only for the 
largest lengths. For a better understanding of this issue, we also attempted to 
perform a more general analysis on the whole considered temperature range 
by introducing a finite correlation length, with the aim of partially taking 
into account finite size effects. Results in the region T ~ T c point towards a 
c < 2 value and are our best numerical evidence of a second order transition. 
Anyway, it should be stressed that we get a quite large c ~ 1.9 and we can 
not rule out the possibility of a sharper (first order) transition. 

During completion of this work, results appeared on the homogeneous 
Poland- Scheraga model [35] that seem to show how this model with the ap- 
propriate parameter values and the lattice SAW-DNA model are equivalent, 
even though the lattice model displays strong finite size corrections to scal- 
ing. Therefore a different kind of analysis on the disordered Poland- Scheraga 
model could help in better understanding the situation. To this extent, the 
recent study by Garel and Monthus [36] points towards the direction of a first 
order transition also in the presence of disorder, though the usual scaling laws 
seem to be not fulfilled. On the other hand, a very recent theoretical work 
[37] predicts that disorder should be relevant in this kind of models and that 
one should find a second order (or smoother) transition. In a nutshell, the 
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situation appears far from being clarified and these simple models for DNA 
denaturation transition seem to deserve a careful analysis from the statistical 
mechanics point of view. 

Finally, it should be pointed out that even if the transition is smoother in 
the presence of disorder the obtained value c ~ 1.9 for the considered inter- 
action energies is definitely higher than the value c = 1.76276 [32] which one 
gets in the homogeneous case when self-avoidance between different loops and 
segments is neglected. In the hypothesis that c depends on the energy values 
one would expect, in the more realistic case of ecc closer to cat, an higher 
value of c possibly indistinguishable from the case of a first order transition 
with c > 2. In this sense the model seems relevant for the experimentally 
observed DNA denaturation, and which value of c one should use in predic- 
tions based on Poland-Scheraga models as realistic as possible seems to be 
an open question. 



Acknowledgments 

It is a pleasure to thank Serena Causo and Peter Grassberger, who were 
involved in the beginning of this study. I also acknowledge stimulating dis- 
cussions with Enrico Carlon, David Mukamel, Henri Orland, Enzo Orlandini, 
Attilio Stella and Edouard Yeramian. I am moreover indebted to Alain Bil- 
loire for carefully reading the manuscript and to Thomas Garel and Cecile 
Monthus for suggestions and discussion of their results on the disordered 
Poland-Scheraga model with c = 2.15 [36J. This work was partly funded by 
a Marie Curie (EC) fellowship (contract HPMF-CT-2001-01504). 

References 

[1] R.M. Wartell and A.S. Benight, Phys. Rep. 126, 67 (1985). 

[2] D. Poland and H.A. Scheraga, J. Chem. Phys. 45, 1456 (1966); 45, 1464 
(1966). For a review see D. Poland and H.A. Scheraga (eds.), Theory of 
Helix-Coil Transitions in Biopolymers, (Academic, New York, 1970). 



23 



[3] M. Peyrard and A.R. Bishop, Phys. Rev. Lett. 62, 2755 (1989); T. Daux- 
ois, M. Peyrard and A.R. Bishop Phys Rev. E 47, 684 (1993). S. Ares 
el al, Phys. Rev. Lett. 94, 035504 (2005). 



[4] 

[5] 

[6] 

[7] 

[8 
[9 

[10 

[11 

[12 
[13 
[14 
[15 

[16 

[17 



S. Cocco and R. Monasson, Phys. Rev. Lett. 83, 5178 (1999). 

V. Ivanov, Y. Zeng and G. Zocchi, Phys. Rev. E 70, 051907 (2004); V. 
Ivanov, D. Piontkovski and G. Zocchi, Phys. Rev. E 71, 041909 (2005). 

N. Theodorakopoulos, T. Dauxois and M. Peyrard, Phys. Rev. Lett. 85, 
6 (2000). 

M.S. Causo, B. Coluzzi, and P. Grassberger, Phys. Rev. E 62, 3958 
(2000). 

Y. Kafri, D. Mukamel, and L. Peliti, Phys. Rev. Lett. 85, 4988 (2000). 

E. Carlon, E. Orlandini, and A.L. Stella, Phys. Rev. Lett. 88, 198101 
(2002). 

T. Garel, C. Monthus, and H. Orland, Europhys. Lett. 55, 132 (2001); 
S.M. Bhattacharjee, Europhys. Lett. 57, 772 (2002); T. Garel, C. Mon- 
thus, and H. Orland, Europhys. Lett. 57, 774 (2002). 

For a recent review in which the effects of self-avoidance in Poland- 
Scheraga models are critically discussed see C. Richard and A. J. 
Guttman, J. Stat. Phys. 115, 943 (2004). 

R. D. Blake et al, Bioinformatics, 15, 370 (1999). 

M.E. Fisher, J. Chem. Phys. 45, 1469 (1966). 

Y. Kafri, D. Mukamel, and L. Peliti, Eur. Phys. J. B 27, 135 (2002). 

B. Duplantier, Phys. Rev. Lett. 57, 941 (1986); J. Stat. Phys. 54, 581 
(1989). 

M. Baiesi, E. Carlon, E. Orlandini, and A. L. Stella, Eur. Phys. J. B 
29, 129 (2002). 

M. Baiesi, E. Carlon, and A. L. Stella, Phys. Rev. £66, 021804 (2002); 



24 



[18] M. Baiesi, E. Carlon, Y. Kafri, D. Mukamel, E. Orlandini, and A. L. 
Stella, Phys. Rev. E67, 021911 (2003). 

[19] E. Orlandini, S. M. Bhattacharjee, D. Marenduzzo, A. Maritan, and F. 
Seno, J. Phys. A 34, L751 (2001). 

[20] R. Blossey and E. Carlon, Phys. Rev. E 68, 061911 (2003). 

[21] A. Hanke and R. Metzler, Phys. Rev. Lett. 90, 159801 (2003); Y. Kafri, 
D. Mukamel, and L. Peliti, Phys. Rev. Lett. 90 159802 (2003). 

[22] T. Garel and H. Orland, Biopolymers 75, 453 (2004). 

[23] D. Cule and T. Hwa, Phys. Rev. Lett. 79, 2375 (1997). 

[24] D.K. Lubensky and D.R. Nelson, Phys. Rev. Lett. 85, 1572 (2000); Phys. 
Rev. £65, 031917 (2002). 

[25] A. B. Harris, J. Phys. C7, 1671 (1974). 

[26] M. Aizenman and J. Wehr, Phys. Rev. Lett. 62, 2503 (1989); Phys. Rev. 
Lett. 64, 1311(E) (1990). 

[27] A. N. Berker, Physica ,4194, 72 (1993). 

[28] E. Eisenriegler, K. Kremer, and K. Binder, J. Chem. Phys. 77, 6296 
(1982). 

[29] P. Grassberger, Phys. Rev. E 56, 3682 (1997). 

[30] H. Frauenkron, M.S. Causo, and P. Grassberger, Phys. Rev. E 59, R16 
(1999). 

[31] S. Caracciolo, M.S. Causo, P. Grassberger, and A. Pelissetto, J. Phys. 
A 32, 2931 (1999). 

[32] P. Belohorec and B. Nickel, Accurate Universal and Two-parameter 
Model Results from a Monte-Carlo Renormalization Group Study, 
Preprint - University of Guelph (1997). 

[33] S. Wiseman and E. Domany, Phys. Rev. E 52, 3469 (1995); Phys. Rev. 
Lett. 81, 22 (1998); Phys. Rev. E 58, 2938 (1998). 



25 



[34] C. Tebaldi, M. De Menech, and A.L. Stella, Phys. Rev. Lett. 83, 3952 
(1999). 

[35] L. Schafer, Can Finite Size Effects in the Poland- Scheraga Model 
Explain Simulations of a Simple Model for DNA Denaturation?, 
cond-mat/0502668 

[36] T. Garel and C. Monthus, J. Stat. Mech. P06004 (2005). 

[37] G. Giacomin and F.L. Toninelli, Smoothening effect of quenched disorder 
on polymer depinning transition, |math.PR/050643lj 



26 



